from util import *
from glob import glob
/usr/local/lib/python3.10/dist-packages/geopandas/_compat.py:123: UserWarning: The Shapely GEOS version (3.11.1-CAPI-1.17.1) is incompatible with the GEOS version PyGEOS was compiled with (3.10.4-CAPI-1.16.2). Conversions between both will be slow. warnings.warn( /home/nyou045/git/CliffErosion_Projections/util.py:3: UserWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas will still use PyGEOS by default for now. To force to use and test Shapely 2.0, you have to set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas: import os os.environ['USE_PYGEOS'] = '0' import geopandas In a future release, GeoPandas will switch to using Shapely by default. If you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html). import geopandas as gpd # Geospatial data
df = load_AOIs()
df
| Taranaki | AOI | SSP 4.5 (p50) | SSP 4.5 (p83) | SSP 8.5 (p50) | SSP 8.5 (p83) | Rate SSP 4.5 (p50) | Rate SSP 4.5 (p83) | Rate SSP 8.5 (p50) | Rate SSP 8.5 (p83) | match | match_score | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 7 | NORTH | TongaporutuRiver | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | TongapurutuRiverCliffs | 93.750000 |
| 11 | SOUTH | HangatahuaRiver_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | HangatahuRiver_South | 97.435897 |
| 21 | SOUTH | Rahotu | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | Rahotu | 100.000000 |
| 20 | SOUTH | Pihama | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Pihama | 100.000000 |
| 19 | SOUTH | OpunakeBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OpunakeBeachCliffs | 100.000000 |
| 18 | SOUTH | OhaweBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | OhaweBeach | 100.000000 |
| 17 | SOUTH | Oeo | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oeo | 100.000000 |
| 16 | SOUTH | Manutahi | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Manutahi | 100.000000 |
| 15 | SOUTH | ManaBay | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | ManaBayCliffs | 100.000000 |
| 14 | SOUTH | KaupokonuiBeach | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | KaupokonuiBeach | 100.000000 |
| 13 | SOUTH | Kakaramea | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Kakaramea | 100.000000 |
| 12 | SOUTH | Hawera_WaihiBeach | 0.57 | 0.78 | 0.83 | 1.10 | 0.0057 | 0.0078 | 0.0083 | 0.0110 | Hawera_WaihiBeach | 100.000000 |
| 0 | NORTH | Mohakatino_River | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | MohakatinoRiver | 100.000000 |
| 10 | SOUTH | CapeEgmont | 0.58 | 0.78 | 0.84 | 1.10 | 0.0058 | 0.0078 | 0.0084 | 0.0110 | CapeEgmont | 100.000000 |
| 9 | NORTH | Waitara | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Waitara | 100.000000 |
| 8 | NORTH | UrenuiRiver_North | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | UrenuiRiverNorthCliffs | 100.000000 |
| 6 | NORTH | PariokariwaPoint | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | PariokariwaPointCliffs | 100.000000 |
| 5 | NORTH | Onaero | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OnaeroCliff | 100.000000 |
| 4 | NORTH | Oakura_South | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | OakuraSouth | 100.000000 |
| 3 | NORTH | Oakura | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | Oakura | 100.000000 |
| 2 | NORTH | NewPlymouth_Airport | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthAirport | 100.000000 |
| 1 | NORTH | NewPlymouth | 0.57 | 0.78 | 0.84 | 1.10 | 0.0057 | 0.0078 | 0.0084 | 0.0110 | NewPlymouthCliffs | 100.000000 |
| 22 | SOUTH | WainuiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WainuiBeach | 100.000000 |
| 23 | SOUTH | WaipipiBeach | 0.57 | 0.78 | 0.83 | 1.09 | 0.0057 | 0.0078 | 0.0083 | 0.0109 | WaipipiBeachCliffs | 100.000000 |
site = df.match.sample(1).iloc[0]
site
'Oakura'
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
if site == "ManaBayCliffs":
print("Flipping")
for k, v in transect_metadata.items():
transect_metadata[k]["Azimuth"] = v["Azimuth"] + 180
results = predict(gdf, transect_metadata)
results
| TransectID | BaselineID | group | Year | ocean_point | linear_model_point | sqrt_model_point | BH_model_point | Sunamura_model_point | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 8 | 1 | 0 | 2100 | POINT (1682577.3607581342 5671120.99636598) | POINT (1682903.9877330572 5670743.528687456) | POINT (1682904.7560665228 5670742.640760029) | POINT (1682880.738653313 5670770.396572942) | POINT (1682903.9877330572 5670743.528687456) |
| 1 | 9 | 1 | 0 | 2100 | POINT (1682590.3312759965 5671130.686919325) | POINT (1682916.521771773 5670744.2823856585) | POINT (1682927.4639775786 5670731.32027307) | POINT (1682893.6026720947 5670771.432300012) | POINT (1682916.521771773 5670744.2823856585) |
| 2 | 10 | 1 | 0 | 2100 | POINT (1682598.6770314383 5671138.73452017) | POINT (1682927.0276047606 5670747.144026966) | POINT (1682945.154076735 5670725.5264186) | POINT (1682904.1986634596 5670774.369794513) | POINT (1682927.0276047606 5670747.144026966) |
| 3 | 11 | 1 | 0 | 2100 | POINT (1682604.6288549996 5671145.294702299) | POINT (1682937.4481723008 5670749.638776907) | POINT (1682964.3080712678 5670717.70770311) | POINT (1682914.576493046 5670776.828651263) | POINT (1682937.4481723008 5670749.638776907) |
| 4 | 12 | 1 | 0 | 2100 | POINT (1682607.6075595547 5671154.595712404) | POINT (1682942.7813369657 5670757.408985613) | POINT (1682972.087082433 5670722.681170065) | POINT (1682919.8669761904 5670784.56289969) | POINT (1682942.7813369657 5670757.408985613) |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 97 | 187 | 4 | 7 | 2100 | POINT (1685286.8791120881 5673060.735761237) | POINT (1685696.9093590495 5672775.544998921) | POINT (1685703.3645548837 5672771.055178392) | POINT (1685667.740752512 5672795.832812311) | POINT (1685696.9093590495 5672775.544998921) |
| 98 | 188 | 4 | 7 | 2100 | POINT (1685286.8521959984 5673068.408532419) | POINT (1685700.9186097626 5672784.572777532) | POINT (1685707.2977494132 5672780.199981628) | POINT (1685671.612585248 5672804.661578058) | POINT (1685700.9186097626 5672784.572777532) |
| 99 | 189 | 4 | 7 | 2100 | POINT (1685287.0361195225 5673075.640930105) | POINT (1685706.5802912419 5672792.548262031) | POINT (1685720.9431551918 5672782.856740254) | POINT (1685677.127796816 5672812.421700373) | POINT (1685706.5802912419 5672792.548262031) |
| 100 | 190 | 4 | 7 | 2100 | POINT (1685289.1618428899 5673080.93252381) | POINT (1685743.416889592 5672780.039583565) | POINT (1685834.9939295463 5672719.380081791) | POINT (1685713.7955143917 5672799.660417153) | POINT (1685743.416889592 5672780.039583565) |
| 101 | 191 | 4 | 7 | 2100 | POINT (1685294.8243254684 5673083.481260104) | POINT (1685719.8676822786 5672807.560334522) | POINT (1685733.991497258 5672798.391726098) | POINT (1685690.0660794459 5672826.906327597) | POINT (1685719.8676822786 5672807.560334522) |
102 rows × 9 columns
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
poly = prediction_results_to_polygon(results)
poly.explore(tiles="Esri.WorldImagery")
site = df.match.sample(1).iloc[0]
site = "CapeEgmont"
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
m = gpd.read_file(f"Shapefiles/{site}_TransectLines.shp").explore("TransectID", tiles="Esri.WorldImagery")
results = predict(gdf, transect_metadata)
results["geometry"] = results["linear_model_point"]
pd.concat([gdf, results])[["Year", "geometry"]].explore("Year", m=m)
gpd.GeoSeries(results.ocean_point, crs=2193).explore(m=m, color="blue")
#gpd.GeoSeries(poly, crs=2193).explore(color="pink", m=m)
results.set_geometry(f"linear_model_point", inplace=True, crs=2193)
polygon = prediction_results_to_polygon(results)
polygon.explore(m=m)
m
CapeEgmont
run_all_parallel()
0%| | 0/24 [00:00<?, ?it/s]
Flipping
samples = pd.concat(gpd.read_file(f) for f in glob("Projected_Shoreline_Polygons/*_linear_Rate SSP 4.5 (p50).shp"))
samples.explore(tiles="Esri.WorldImagery")
# Compare different SLR rates
site = df.match.sample(1).iloc[0]
print(site)
gdf = gpd.read_file(f"Shapefiles/{site}_intersects.shp")
gdf.crs = 2193
gdf = enrich_df(gdf)
transect_metadata = get_transect_metadata(f"Shapefiles/{site}_TransectLines.shp")
low_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.007)
high_proj_slr = predict(gdf, transect_metadata, Proj_SLR=0.015)
m = gpd.GeoSeries(low_proj_slr.sqrt_model_point, crs=2193).explore(color="blue", tiles="Esri.WorldImagery")
gpd.GeoSeries(high_proj_slr.sqrt_model_point, crs=2193).explore(color="green", m=m)
WaipipiBeachCliffs